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Abstract 

In this paper, we propose designing transportation network topology and traffic distribution under fluctuating conditions 
using a bio-inspired algorithm. The algorithm is inspired by the adaptive behavior observed in an amoeba-like organism, 
plasmodial slime mold, more formally known as Plasmodium of Physarum plycephalum. This organism forms a 
transportation network to distribute its protoplasm, the fluidic contents of its cell, throughout its large cell body. In this 
process, the diameter of the transportation tubes adapts to the flux of the protoplasm. The Physarum algorithm, which 
mimics this adaptive behavior, has been widely applied to complex problems, such as maze solving and designing the 
topology of railroad grids, under static conditions. However, in most situations, environmental conditions fluctuate; for 
example, in power grids, the consumption of electric power shows daily, weekly, and annual periodicity depending on the 
lifestyles or the business needs of the individual consumers. This paper studies the design of network topology and traffic 
distribution with oscillatory input and output traffic flows. The network topology proposed by the Physarum algorithm is 
controlled by a parameter of the adaptation process of the tubes. We observe various rich topologies such as complete 
mesh, partial mesh, Y-shaped, and V-shaped networks depending on this adaptation parameter and evaluate them on the 
basis of three performance functions: loss, cost, and vulnerability. Our results indicate that consideration of the oscillatory 
conditions and the phase-lags in the multiple outputs of the network is important: The building and/or maintenance cost of 
the network can be reduced by introducing the oscillating condition, and when the phase-lag among the outputs is large, 
the transportation loss can also be reduced. We use stability analysis to reveal how the system exhibits various topologies 
depending on the parameter. 
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Introduction 

Transportation networks such as power grids are, in general, 
designed under certain static supply-demand conditions. However, 
in most situations, whether the network is that of nature or a man- 
made system, the inputs/ outputs into/ from the networks fluctuate 
rather than remain constantly static. One such example is an ant 
foraging trail network, in which ants cannot constantly prey upon 
their foods because the activities of the prey animals fluctuate daily 
or seasonally. The feature also holds for man-made networks. For 
instance, the number of passengers commuting by rail is 
maximized in the mornings and evenings and the peak times shift 
among stations in suburbs and city areas on weekdays. Addition- 
ally, the transportation patterns on weekends are quite different 
from those on weekdays. The second man-made example is power 
grids. The pattern of electricity consumption is distributed 
according to the lifestyles or business style of consumers, which 
was recendy confirmed using clustering analysis on a town in 
Japan [1]. More specifically, the consumption pattern fluctuates 
daily, weekly, and seasonally, and the peak time depends on the 
consumers. 

Optimization of networks under fluctuating conditions is difficult 
to be conducted in a straightforward manner by conventional 



methods within linear- and nonlinear-programming frameworks. In 
this paper, we propose designing traffic distribution in networks 
under fluctuating conditions using an algorithm inspired by the 
organism Physarum. 

The Physarum algorithm, which mimics the shortest path-finding 
behavior of the plasmodial slime mold organism [2], formally 
called Physarum polycephalum, was developed by Tero et al. [3] . The 
Plasmodium of Physarum is a giant amoeba-like multinucleated 
unicellular organism. It contains thousands of nuclei, so the cell 
size can get very large, ranging from 10 jim to 1 m. To distribute 
protoplasm, including nutrients, oxygen, and organelles, through- 
out this large cell body, the organism has developed a peculiar 
transportation network consisting of tubular structures. The 
diameter of the tubes adapts to the flux of the protoplasm: The 
tubes on the paths connecting multiple food sites become thick in 
accordance with the growth of the protoplasmic flow, while the 
other paths become thin and finally disappear when there is little 
or no flow. Consequendy, the organism is able to generate the 
shortest paths connecting multiple food sites [2,4]. The Physarum 
algorithm, which mimics the adaptive behavior of the tubes, has 
been widely applied to complex problems such as maze solving 
[2], design of the topology and transportation distribution of 
railroad grids [5,6] and highway networks [7], and path formation 
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in wireless sensor networks [8]. Although, in the above examples, it 
was applied under static conditions, the algorithm can also be 
applied under fluctuating conditions owing to its adaptive 
behavior. 

This paper studies the design of network topology and traffic 
distribution under oscillating conditions, which is the simplest type 
of fluctuating environment. The network consists of nodes and 
links, which, in power grids for example, correspond to consumers, 
power plants, electric poles, and power lines. A multiplicity of 
consumers uses electricity with daily periodicity (oscillating 
condition). The peak consumption times vary according to the 
consumers, and are defined by phase lags. 

In the Methods section, we outiine how the Phy.sarum algorithm 
is modified to deal with problems involving oscillating conditions 
and define performance functions. We then present the network 
designs under oscillating conditions proposed using the Physarum 
algorithm and evaluate them using our performance functions, in 
the Results section. In the Discussion section, a stability analysis 
for a simple network is considered in our discussion of the numeric 
result. Finally, we discuss the effect of the oscillating condition and 
the phase lags. 

Methods 

Physarum Algorithm 

In this section, we modify the original Phy.sarum algorithm [3] to 
deal with the example network depicted in Fig. 1 . The shaded and 
unshaded large circles, respectively, represent nodes for input 
(denoted as in) and output (denoted as OUt\y) of transported 
materials, such as protoplasm in Physarum, current in power grids, 
and people in railroad grids. The link Ig connecting the nodes i 
and j has the following properties: length Ly, conductivity Dg, and 
traffic volume flux Qg, Their meanings in each application are 
summarized in Table 1. 



Table 1. Correspondence of variables. 





Variables 


General 


Physarum 


Power Grid 


Ui 


length of link 


length of tube 


length of electric wire 




conductivity 


tube thickness 

(x4> 


electric conductivity 


Qv 


flux of traffic 
volume 


flux of protoplasm 


current 



The radius of tubes and wires are represented with r$ 
doi:1 0.1 371 /journal.pone.0089231 .tOOl 

The flux at each node conserves 



£2,7 = 



o. 



lin, 



j^in, outi^, 
j = in, 

j = OUt\2, 



(1) 



where /,„ and I ou t l 2 are the fluxes at in and out\ 2, respectively. 
The total flux to/from the system should be balanced: 



I in Iout\ ~t~ Iout2 • 



The flux Qij is given by 



A, 

Qy= ,'iPi-Pj) 

J-'ii 



(2) 



(3) 



where />,- and pj represent, respectively, the pressure at nodes i and 
j. Substituting Eq. (3) for Eq. (1), />, is obtained under the given Ly 



outi i 




Figure 1. Network topology used by the Physarum algorithm for numerical calculation. 

doi:1 0.1 371 /journal.pone.0089231 .g001 
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and Dy. Qy is then calculated using Eq. (3) again. In the numerical 
calculations, Ly = 1.0 is set at all links. 

As mentioned in the introductory section, the conductance Dy 
adapts to flux. Therefore, the conductance is assumed to evolve 
according to the following differential equation: 

d ^=f(\Q i j\)-D ij , (4) 

meaning that the tube grows depending on the flux (the first term 
on the right hand side of the equation) while it degenerates (the 
second term). It is natural in biological systems for the growth rate 
to be saturated by an upper limit. Thus, function /(•) can be 
defined as a sigmoid function: 



/fl&D-T^ji • ( 5 ) 

This function is widely found in biological cooperative processes 
[9] . The parameter fi is the key parameter governing the dynamics 
of this system. When [i> 1, the tube grows only slighdy when the 
flow is extremely weak, although the growth speed is accelerated 
when it once starts to grow, then it is finally saturated at one. 
When [1=1, the function is categorized into the Michaelis- 
Menten type, which represents the simplest enzymatic reaction. 
When ji< l,/( ) represents fast initial growth and slow saturation, 
suggesting no meaning related to biological processes. 

Application to Networks with Oscillating Conditions 

The outputs at out\^ are assumed to oscillate as follows: 

I 0 ut x = 1 + sin wt, I ou , 2 = 1 + sin (wt+(f)) , (6) 

where O) is angular frequency and <f> represents a phase lag 
between the outputs. In Eq. (6), w = 2nx 10 2 is assumed so that 
the period of oscillation T = 2n/ai= 10~ 2 is small enough to the 
time constant of the degeneration process of Dy, which is 
estimated as one. We confirmed that the outline of the results in 
this paper is valid over the frequency range 
27ixlO~ 1 <(D<27rxl0 7 , namely, over the period range 
1(T 7 < T< 10 1 (see §1 in File SI for details). 

Converged Value of Dy 

We repeat the computation of Eq. (4) with Eqs. (l)-(3), and Eq. 
(6) until Dy converges within a certain accuracy. In fact, Dy 
continues to oscillate slighdy even after a long evolution period 
(Figs. A and B in §1 in File SI). Therefore, the completion of the 
convergence is judged according to the following criterion. 

First, a variation amount of Dy at time t is defined as follows: 

e(t)=^{D iJ (t+nT)-Dy(t)} 2 , (7) 

v 

where «(=100) is the number of cycles in the input/ output 
oscillation, and m is the total number of links. The convergence of 
Dy is judged to be complete when e(i) < e c = 1 0 ~ 26 (the value of e c 
is set by the reasoning below). Consequendy, the averaged value 
Dy at time t + nT is calculated using the following definition: 



Dy(t)=^ T Dy( S )ds. (8) 

After the convergence is ascertained, Dy(t) is denoted as Dy. The 
link is removed (Dy set to zero) when Dy becomes less than a 
certain threshold -D m ; n = 10~ 13 . The value is sufficiently smaller 
than the order of the maximum value of Dy (~10 _1 ). Finally, the 
network topologies and traffic distributions (magnitude of 
conductances) recommended by the Physarum algorithm are 
obtained (Figs. 2 and 3). Note that the threshold value for 
judgment of Z)^ -convergence e c is sufficiendy smaller than D m [ n . 

Performance Functions 

We now introduce three performance functions to evaluate the 
performance of the networks recommended by the Physarum 
algorithm: power or transportation loss, building and/or mainte- 
nance cost, and vulnerability in network topology. 

Loss P is defined using an analogy to electric energy loss, which is 
calculated with (current) 2 x(electric resistance) in a wire. Conseq- 
uently, the loss for alink/y is defined as Qy multiplied by Ly/Dy 
(see also §2 in File SI). The total loss for the network is calculated 
by summing thelossforeachlinkoverallthelinksasfollows: 




where the loss is averaged over a period of input/output oscillation 
because Qy oscillates. 

Cost B is that for building and/or maintaining a network, which 
is expected to be proportional to the total volume of the network. 
Because the cross section of each link is proportional to Dy in the 
case of power grids, as described in Table 1, B is defined as 
follows: 

B >>A- (10) 

ij 

Note that B= Y^y Ly\J~D~y should be adopted when consider- 
ing the original Physarum network because the relation between 
conductivity and a tube of radius r-y is described as Dyazry (see 
also Table 1) [3,6]. 

Vulnerability V is defined as the probability that the connection 
out\ or out2 from in is divided when one of the links in the network 
is randomly deleted. The deletion frequency is assumed to be 
proportional to the length of the link when Ly is not homoge- 
neous, where the probability is normalized by the total length of 
the network, L tot . Consequently, the vulnerability is defined as 
follows: 

v =J2 G v J r L ' (n) 

y ^tOt 

where disconnectivity Gy for a link ly is set to one if transportation 
flows out from in can reach neither out\ nor outi\ otherwise, it is 
set to zero. 
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Figure 2. Dependence of network topology under constant input and output. Initial values of D u were either set as homogeneous 
(Ay = 10 for all links) or were distributed according to a normal distribution with mean 1.0 and standard deviation 0.1 . Solid and dashed lines of the 
network diagrams denote surviving and removed links, respectively. A Complete mesh (type 1), 0</j < 1. B Partial mesh (types 2-5), 1 </x< 1.4. C V- 
shaped network (types 6 and 7), fi > 1.4. D Y-shaped network (type 8), /i > 1.8. When the initial conditions of Dg are exactly homogeneous, the V- 
shaped network appears in the range of 1.4 </,i<2.2 and the Y-shaped network appears in the range of /i> 2.3. Type numbers correspond to those 
of Figs. 3 and SI. 

doi:1 0.1 371 /journal.pone.0089231 .g002 



Results 

We considered two cases, constant and oscillating flux at input 
and output nodes, and evaluated the network topologies and traffic 
distributions recommended by the Physarum algorithm using the 
three performance functions. 

Constant Condition 

Before capturing the effect of oscillatory input/output on the 
network design, we tested the effect with constant input/output. 
We set the fluxes to constant values, /,„ =2, I ou t x ,ouh = 1, in Eq. (1). 
The numerical calculation started from a homogeneous initial 
condition of Dy = 1 .0 or a non-homogeneous condition according 
to normal distribution with mean 1.0 and standard deviation 0.1. 
We observed eight types of network topologies in the parameter 
range Q<fi<5.0 as shown in Fig. 2. 





o 
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nu 
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CD 
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Figure 3. Network types calculated with oscillating inputs and 
outputs. A type number is assigned to each topology. The full list of 
network topologies is represented in Fig. S1 . The data are those for the 
homogeneous initial conditions of Dy. The plots of mesh, partial mesh, 
V-shaped and Y-shaped networks, are colored in black, red, blue, and 
green, respectively. The dependence of type of partial mesh on ^ is 
shown in Fig. S2. 

doi:1 0.1 371 /journal.pone.0089231 .g003 



The network topology changes from dense to sparse depending 
on fi. When fi is smaller [fi < 1.4), the network forms a mesh 
accompanied by circular structures (Figs. 2A and 2B). When \i 
becomes larger (jj. > 1 .4), the network forms a tree structure 
(Figs. 2C and 2D). The mesh networks are categorized into two 
types, complete mesh (type 1; Fig. 2A) and partial mesh (types 2-5; 
Fig. 2B). The tree networks are categorized into two types, V- 
shaped (types 6 and 7; Fig. 2C) and Y-shaped (type 8; Fig. 2D) 
networks. The Y-shaped networks appear when ft > 1 .8. The paths 
from the input are partially shared in the Y-shaped network, while 
they are directly connected to the two outputs in the V-shaped 
network. 

Oscillating Condition 

We set the input/output flux oscillating using the definition in 
Eq. (6). The numerical calculation started from a homogeneous 
initial condition of Dy = 1.0 or non-homogeneous conditions 
according to normal distribution with mean 1.0 and standard 
deviation 0.1. We observed 20 types of network topologies in the 
parameter range 0<^i<5.0, as shown in Figs. 3, SI and S2. In 
this case, the dependence of network topology on fi is similar to 
that of the constant condition: when 0 < fi < 1 , complete mesh 
(type 1) appeared. As fi increased over 1, the topology changes to 
partial mesh (types 2-5 and 9-19). Finally, when /J > 1.5, V- 
shaped (types 6 and 7) or Y-shaped (type 8) networks were 
observed. It should be noted that the variation of the topologies 
becomes broader than in the case of the constant condition when 
1 S 1-8: a variety of partial meshes, i.e., networks of types 9-19 
besides types 2-5, were observed. 

The network topology depends not only on fi but also on phase 
lag (f) and on the initial conditions of Dy. The characteristics are 
particularly evident in fi > 1 . Figure 4 shows the network types 
observed according to fi, <f>, and the initial conditions of Dy. For 
fi < 1.7, primarily partial meshes were observed (see Fig. S2 for 
details). Dependence of the topology on 0 can be seen more clearly 
when fi > 1.7: the V-shaped network is more frequently observed 
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Figure 4. Relation between network types and the parameters ii and (fi under oscillating conditions. A Homogeneous initial condition of 
Djj—l.Q. B D Examples of non-homogeneous initial condition of Dy-. Initial values of £>,y were distributed according to a normal distribution with 
mean 1.0 and standard deviation 0.1. Black, gray, and white squares denote partial mesh, V-shaped and Y-shaped networks, respectively. The specific 
type-number of partial mesh depends on both parameters /( and (p (Fig. S2), and also on the initial condition of £>,y- which is not shown here in detail. 
doi:10.1371/journal.pone.0089231.g004 
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Figure 5. Performance depending on parameters /i and 0. A Loss P. B Cost B. C Vulnerability V. Circles, triangles and squares denote 
performances when (j) = 0,Ti/2,n, respectively. The crosses in C denote the performances of the constant condition. The data are those for the 
homogeneous initial conditions of Dg. The case starting from non-homogeneous initial conditions is demonstrated in §3 in File S1. 
doi:1 0.1 371 /journal.pone.0089231 .g005 
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Figure 6. Comparison of the performances for the networks designed under constant and oscillatory conditions. A Ratio in loss, 
P 0 /Pc- B Ratio in cost, B 0 /B c . Circles, triangles, and squares respectively denote (h=Q,%/2,n. The data are those for the homogeneous initial 
conditions of £>,y. 

doi:1 0.1 371 /journal.pone.0089231 .g006 



when the two outputs are in phase (<fi~0) and the observation 
ratio of the Y-shaped network increases accordingly as the lag 
approaches anti-phase Note that the dependence of the 

topology on <j) is also subject to the initial distribution of Dy in 
detail (compare the diagrams A of homogeneous condition and B 
D of three different non-homogeneous conditions in Fig. 4). 



Evaluation of the Networks 

Figure 5 shows the performances P, B, and V estimated for 
each combination of parameters fi and <f), where each network is 
calculated from the homogeneous initial conditions of Z),/. Smaller 
values mean better performances in these analyses. Loss P 
increases until around j.i = 1.7, then slightly decreases, irrespective 
of (f>, as shown in Fig. 5A. Notably, P for <f> = n is clearly always 
smaller than those for (p = 0 and n/2. The discontinuity in the plots 
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Figure 7. The simple network used for stability analysis. The link 

lengths were set as L\=L2=L^ = l.O and Li, = L^ =2.0 so that any 
path length from in to out\j is 2. 
doi:1 0.1 371 /journal.pone.0089231 .g007 

for (f) = 0 when 4</i<4.4 is caused by the discontinuous change 
of network topology. Cost B decreases rapidly until around /i= 1, 
then it becomes almost constant, as shown in Fig. 5B. Vulnera- 
bility V equals 0 when fi< 1.5, as shown in Fig. 5C because the 
network includes circular structures (Figs. 2A and B). As fi exceeds 
around 1.5, V jumps to 1.0 because the network includes no 
circular structure. In conclusion, the network is well balanced at 
1 S 5 1-5. The results for the non-homogeneous initial conditions 
of Dy are valid for virtually the same feature as in the case for the 
homogeneous conditions (see §3 in File SI for details). 

Benefit Derived from the Introduction of Oscillatory 
Condition 

To investigate the benefit derived from the introduction of the 
oscillatory condition, we calculated the ratio of the performances 




Figure 8. Twelve equilibria for the network in Fig. 7 represented i 

shaped network, H-J V-shaped network, K and L the others. 
doi:1 0.1 371 /journal.pone.0089231 .g008 



between the constant and oscillatory input/output, as shown in 
Fig. 6. Note that the performance P c and B c were estimated with 
oscillatory input/output against the networks obtained under 
constant condition by the Phy.samm algorithm. The performances 
P 0 and B 0 were estimated with oscillatory input/output against 
the networks obtained under the oscillatory condition, which are 
the same as those of Fig. 5. 

A ratio with value smaller than 1.0 suggests that the 
performance of the network considering the oscillatory condition 
is better. Loss P for the oscillatory condition is better than that for 
the constant condition only when 1.4 < jU < 1.5. In contrast, cost B 
almost always shows better performance in the oscillatory 
condition. The cost can be reduced to about 80% in the best 
performance. The effect of vulnerability is captured in Fig. 5C: 
Vulnerability is improved by considering the oscillatory condition 
when 1.3 </j< 1.8. 

Discussion 

Stability Analysis of Network Topology 

To understand the parameter dependence of the network 
topology, we conducted stability analyses of network topologies 
and estimation of their basin size against a network with small 
compositions of nodes and links (Fig. 7). 

In this subsection, the notation of link /« is redefined as l/ c . In 
accordance with this definition, the equations for conductances 
D = (Di, ■ ■ ■ ,D$) are rewritten instead of using Eq.(4) as follows: 

^=f(\Qk\)~D k , fe=l,...,5. (12) 




network-topology form. A Complete mesh, B-F partial mesh, G Y- 
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Figure 9. Observation ratio for small network. The network in Fig. 
7 was used. Y, V, and V respectively denote the network topologies 
represented in Figs. 8G : 8H, and 81— J. The parameter /< = 2.0 was set. In 
each calculation against d>, all combinations among = 0.1,0.2, ■ ■ ■ ,1.0 
for all links i= 1,2, ■ ■ ■ ,5 (specifically, a total of 10 5 combinations) were 
tested as initial conditions. The networks such as the ones shown in Fig. 
8K and L were also observed but the observation ratios were extremely 
small, e.g., 0 when Q<(j)<0.1n, 0.002 when (j) = Q.^n, 0.018 when 
(f> = 0.9n, and 0.03 when (j) = n. 
doi:1 0.1 371 /journal.pone.0089231 .g009 



Iout 2 =0.5{1 + sin(oj? + ri>)}, where the magnitude of the input/ 
output flux is set as half of those in Eq. (6) because the network size 
is now reduced. 

In Eq. (12), Dk has two time scales: slow and fast. The fast time 
scale is caused by ft), which gives fluctuations with small amplitude 



to D/c. Accumulation of the small asymmetric fluctuations finally 
derives a slow drift in Dk. The final network topology must be 
determined mainly by the slow dynamics. Therefore, Dk can be 
averaged over a period of the fast dynamics when we focus only on 
slow dynamics, which is denoted as Dk hereafter. The slow 
dynamics of Dk can be written as follows: 



where f k 



dDk -f T> 
Tj 0 



=gk 



f{\Qk(t)\}dt 



(13) 



The steady state of Eq. (13), 



dt 



= 0, is considered then the 



solutions of the equation g£ = 0, namely equilibria, are denoted as 
D* = {D\, ■ ■ ■ ,.Dg) . The magnitudes of individual elements of D* 
determine the topology of the network. Note that Qk, and also fk, 
are a function of Dk owing to Eq.(3). Therefore, we solved 
equation gk = 0 using Newton's method, where fj is obtained by 
numerical integration of /(•) according to the above definition 
using Eqs. (l)-(3), (5). The integration otf{\Qk(t)} over the period 
of output oscillation in Eq. (13) depends on (f> because of the 
nonlinearity of the function (see §4 in File SI for details). 

We obtained 12 equilibria of D*, as summarized in Fig. 8, 
where the topologies are drawn based on the magnitude of the 
elements' values, D\,---,Dy The topologies can be roughly 
classified into complete mesh (Fig. 8A), partial mesh (Figs. 8B— F), 
Y-shaped (Fig. 8G), V-shaped (Figs. 8H^J) networks, and others 
(Figs. 8K and L) similar to those of Fig. 2 and SI. The V-shaped 
network is, furthermore, divided into subcategories: symmetric 
(Fig. 8H, denoted by the V-shaped network in Fig. 9), and 
asymmetric (Figs. 8I^J, denoted by the V'-shaped network in 
Fig. 9). 

We conducted linear stability analysis for each equilibrium D*. 
The Jacobian matrix J of g = (gi, • • • ,gs) is defined using Eq. (13) 
as follows: 



O 



X 
CO 

E 



o 

o 
o 

o 



o Complete mesh (X max I 
A Partial mesh Umax) 
□ Y shape (^ max ) 
x V shape (Xtc) 



0.0 



0.5 



1.0 



• B • • • M- ■ 



1.5 



2.0 



Figure 10. Maximum eigenvalues depending on /xwhen</>=0. Circles, triangles, and squares, respectively, denote /„„„ at the equilibria of 
complete mesh (Fig. 8A), partial mesh (Fig. 8E), and Y-shaped (Fig. 8G). Crosses represent X tc for V-shaped (Fig. 8H) networks. 
doi:1 0.1 371 /journal.pone.0089231. g010 



PLOS ONE | www.plosone.org 



9 



February 2014 | Volume 9 | Issue 2 | e89231 



Network Design with Fluctuating Environment 



Table 2 



. Evaluation of the network type for each item. 



V- 


Network\ 
Evaluation 


Loss Cost 


Vulnerability 


Cascading 


<1.0 


Complete mesh 


A + C 


A + 


C 


1.0-1.4 


Partial mesh 


B A 


A + 


C 


>1.4 


V-shaped or 
Y-shaped 


B or C A + 


C 


A + 



A + : best, A: good, B: acceptable, C: bad. 
doi:1 0.1 371/journal.pone.0089231 .t002 



Jkl 



Sgk 
3D,' 

eh 

3D, 

'§k_ 
3D, 



-.5 



=< 



eh 

8Di 



k,l=\, 

dDk 
3D, 

-1 (*=/), 



(14) 



Because it is difficult to calculate Eq. (14) directly, we estimated 
the Jacobian matrix at D* (denoted as J* hereafter) using the 
following approximate form: 



dfk 
3D, 



D = D* 



/t(D*+AP/)- 
\AD,\ 



-fk(P*) 



(15) 



where AD/ is a vector with an l-th element valued S and the others 
zero, e.g., AD2 = (0,<5,0,0,0). For the numerical calculation, 
r)=10~ 5 was used. We then calculated the eigenvalues for J*, 

-<0, the equi- 



A\, X$. When max[Re(Ai), • • • ,Re(A 5 )] ; 
librium D* is determined as stable. 

The above method is not appropriate to examine whether the V- 
shaped network (Fig. 8H) is globally stable because changing the V- 
shaped network (Fig. 8G) to other network types, such as complete 
or partial mesh (Fig. 8A, D, E, or F), requires at least two additional 
links. In Eq. (14), only a single additional link can be considered. 
Therefore, instead of calculating eigenvalues, we estimate a time 
constant X lc converging to D* from a vicinity. We tested four 
combinations of deviations from the V-shaped equilibrium, 
(5,5,5,5,5), (3,3,6,8,-3), {6,3,5,-6,5), (5,6,5,-6,-6). Finally, 
we defined the maximum time constant as kt c - 

Figure 10 summarizes the dependence of the maximum 
eigenvalues l„ mx (or X tc for the V-shaped network) on the 
parameter \i. The single stable equilibrium, complete mesh, is 
found in the region of 1.0. The complete mesh remains stable 
over fi> 1.0 followed by participation of the Y-shaped, V-shaped, 
and partial mesh networks. The complete and partial meshes 
become unstable when ji exceeds 1.3. The stability change from 
complete mesh, via partial mesh, to Y-shaped or V-shaped 
network resembles that of the larger network (Fig. 3). However, no 
significant difference can be found in the features of the stability 
among different phase-lags 0 = 0, n/2 and n (Fig. S3) while 
appearance of Y-shaped or V-shaped network apparently depends 
on cj) in the larger network, as seen in Fig. 4. The dependence 
would be caused by the difference in the basin sizes between the Y- 
shaped and V-shaped networks. Figure 9 shows the observation 



ratio of the Y-shaped and the V-shaped networks. Both types are 
always observed but the ratio of the Y-shaped network increases in 
accordance with (j). The change in basin size depending on (j) could 
explain the observation that the Y-shaped network is more 
frequently observed in anti-phase lag in the larger network, as seen 
in Fig. 4. 

Summary and Conclusion 

In this paper, we proposed using the Phjsarum algorithm to 
design transportation network topologies and traffic distribution 
under oscillating conditions. The results of numerical experiments 
indicate that this approach is valid and has the following benefits: 

1) Only one parameter \l can control the morphology of the 
network. The client using the network can choose a particular 
parameter according to which they consider to be the most 
important among loss, cost, and vulnerability. 

2) By introducing oscillating condition, building and/ or main- 
tenance cost is reduced to a maximum of 80% that of cases in 
which conditions are static. 

3) Phase lag among outputs results in a wide variety of network 
morphology when f,i> 1 (sigmoidal growth in the conduc- 
tance). 

Table 2 summarizes the first item. Partial mesh can be 
recommended when the client requests a system with loss, cost, 
and vulnerability well-balanced. The third index, vulnerability, 
should be noted when considering power grids. The meshed 
network has a low vulnerability index but it includes loop 
connections, which are prone to cascading failure problems. 
When some nodes or links in a meshed network are damaged, the 
current that would normally go through those links must be 
distributed to the surrounding links. However, if the current goes 
beyond the capacity of the surrounding links, the damage 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 
Phase lag (j) 

Figure 11. Relation between loss P and (f) of small network. 

Circles, triangles, and crosses respectively denote Y-shaped network 
(Fig. 8G), V-shaped network (Fig. 8H), and complete mesh (Fig. 8A). The 
total volume ( = cost) for each network is normalized by that for the 
complete mesh so that the three networks are made with the same 
cost. 

doi:10.1 371/journal.pone.0089231 .g01 1 
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propagates rapidly to the outer surrounding links. This results in 
large-scale blackouts [10,11]. Considering these phenomena, V- 
shaped and Y-shaped networks are recommended rather than 
partial mesh. For railroads and highways, in which cascades need 
not be considered, partial mesh can be recommended. The 
cascading problem was not treated as a performance function in 
this paper because, for the sake of simplicity, the capacity of the 
current for each link was not considered. This will be dealt with in 
future work. 

For the second item, if a client considers the reduction of power 
loss more important than building and maintenance cost, a 
network that is designed under static conditions is recommended. 
The recommendation can be reversed by considering the third 
item, phase lag. Then, the problem of loss can be overcome. 

For the third item, the Y-shaped network is observable more 
frequently than the V-shaped network as the phase lag gets larger 
when fi> 1.4. This topological selection delivers a maximum of 
20% loss reduction to the system. Notably, the loss decreases when 
the lag approaches anti-phase away from in-phase, as shown in 
Fig. 1 1 . This result theoretically supports a justification of the 
"peak shift" action developed in Japan for reducing electric power 
after the Fukushima nuclear disaster in 20 1 1 . The peak shift action 
shifts usage of electricity from on-peak to off-peak periods. This 
allows the electric power consumption in power grid systems to 
flatten during the day and to be reduced in peak periods. By 
introducing this action, the number of standby power plants can 
be reduced: Such the plants, e.g., thermal ones, are on standby to 
regulate power generation flexibly and to avert power shortages in 
peak periods. Our results suggest that peak shift action contributes 
to a reduction in not only the number of standby power plants but 
also in power loss in the grid. 

Natural systems may gain advantages by self-organizing their 
network. Argentine ants are known to make supercolonies, which 
consist of multiple colonies with a single family. They form V- 
shaped or Y-shaped trails connecting the multiple colonies [12]. 
Army ants build dendritic trails-large-scaled Y-shaped branching 
structures [13]. Tao et al. showed, by a computer simulation, that 
virtual ants building Y-shaped trails can gain more food than those 
building V-shaped trails when the foods appear in anti-phase at 
two sites [14]. By considering the number of ants as cost, the result 
can be interpreted as the ants selectively building Y-shaped 
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